# Packages ----------------------------------------------------------------# using a reproducible environmentrenv::restore()# pacman allows to check/install/load packages with a single call# if (!require("pacman")) install.packages("pacman") # already in renv.locklibrary("pacman")# packages to load (and install if needed) -------------------------------pacman::p_load( here, # easy file paths see, # theme_modern and okabeito palette report, # reporting various info labelled, # labelled data# ---- packages specific to this project ---- scales, lme4, car, simr, readxl, openxlsx, statmod, qqplotr, emmeans, latex2exp, ggbeeswarm, ggpubr, patchwork, quarto, easystats,# --- tidyverse # modern R ecosystem)# Custom functions shared across scripts ----------------------------------source(here("scripts/_functions.R"))# Global cosmetic theme ---------------------------------------------------theme_set(theme_modern(base_size =14)) # from see in easystats# setting my favourite palettes as ggplot2 defaultsoptions( ggplot2.discrete.colour = scale_colour_okabeito,ggplot2.discrete.fill = scale_fill_okabeito,ggplot2.continuous.colour = scale_colour_viridis_c,ggplot2.continuous.fill = scale_fill_viridis_c)# Fixing a seed for reproducibility ---------------------------------------set.seed(14051998)# Adding all packages' citations to a .bib --------------------------------knitr::write_bib(c(.packages()), file =here("bibliography/packages.bib"))
1 Rationale
We hypothesised that difference in reaction times (RT) between congruent and incongruent trials (the congruence effect) in the tasks would be modulated by the visual imagery ability of participants, and therefore that we would observe between-group differences in the magnitude of this effect. In a statistical model with categorical predictors, evidence in favour of our hypothesis would be provided by a significant two-way interaction between Group and Congruence condition. Estimating the statistical power to capture such an effect, given an experimental design, a sample size and an effect size, requires simulations from a model close to the one that will be used for the analysis of real data, as well as assumptions about the structure of data. The power estimation procedure requires 3 steps:
Specifying the structure and coefficients of a generative model
Generating a synthetic dataset using the generative model and an expected effect size
Fitting the same model on the simulated data and test the statistical significance of the interaction of interest
Steps 2 and 3 are then repeated as many times as necessary to reach a desired precision of estimation. The proportion of synthetic datasets for which the effect of interest is significant is an estimation of power. Steps 2 and 3 can also be repeated for a range of effect sizes, if the expected effect size is not known exactly which is what we have done here. We detail and justify below every assumption made for the generative model.
2 Generative model
In real datasets, the distributions of RTs are right-skewed and can be modelled with transformations using Gamma or inverse Gaussian distributions. However, such a modelling choice is only an approximation as RTs do not follow exactly any distribution of the exponential family for which generalized models can be fitted. It is difficult to evaluate the impact of such approximation on type I and type II errors. Thus, for power calculation, we generated normally distributed RT and modelled them with a regular mixed model, which simplifies the model parametrization as well as alleviates the computational burden.
The generative model has been created using the simr R package.
Creating a dataset following the structure of our data for the generative model
# First creating a dataframe with no response data, but whose structure reflects the dataset that will be collected. # vectors for subjectsid <-seq(1:150)groups <-rep(c("aphantasia", "control"), 75)# conditionscongruences <-list("congruent", "incongruent")colors <-list("colored", "uncolored")# repeated 4 times eachcongruence =rep(congruence, 4)color =rep(colors, 4)# creating the simulation dfdf <-tibble(id = id, group = groups,congruence =list(congruence),color =list(color) ) |># crossing conditions for 64 trials/subject totalunnest_longer(congruence) |>unnest_longer(color)
Fixed effects: We included the Group (aphantasic, control), Congruence condition (congruent or incongruent) and Color condition (color or uncolored) along with all their two and three way interactions as fixed categorical predictors. Regression coefficients were set in such a way as to produce a pattern of means roughly following our hypotheses, where the control group would be faster solely in the congruent condition.
simr requires regression coefficients in order to simulate data. The functions below transform a pattern of cell means to regressions coefficients of a linear model, using a simple trick: modelling a dataframe of hypothesized cell means, as if each of them was one noise-less observation of a combination of experimental factors, provides the desired coefficient estimates.
Creating regression coefficients through models
# creating dataframes of regression coefficientsextract_coefficients <-function(effectsize) { df_temp <-tribble(~group, ~congruence, ~color, ~y,"aphantasia", "congruent", "colored", -30,"aphantasia", "congruent", "uncolored", 0,"aphantasia", "incongruent", "colored", -30, "aphantasia", "incongruent", "uncolored", 0,"control", "congruent", "colored", -30- effectsize,"control", "congruent", "uncolored", 0- effectsize,"control", "incongruent", "colored", -30,"control", "incongruent", "uncolored", 0 ) lm_coef <-lm(y ~ (group + congruence + color)^3, data = df_temp)# Set numerical errors to 0coef(lm_coef)[abs(coef(lm_coef)) <10^(-10)] <-0# converting to ms upon returnreturn(coef(lm_coef)/1000)}
Random intercepts: We included a random intercept for participants with a standard deviation of 150ms. This number was derived from the literature (see Fucci et al., 2023) - however, previous simulation studies have shown that the magnitude of random intercepts has no effect on statistical power whatsoever (Brysbaert & Stevens, n.d.), a result that we could verify in our own simulations.
Random slopes: We included random by-participant slopes for Congruence and Color, thus specifying the full model given our factors and resulting in the most conservative power estimates. Standard deviations of random slopes are virtually never reported in scientific publications; we set them at 50ms, which is in the same range as the maximum expected population-average effect size.
Correlations between random effects were all set to 0, but an alternative value (0.5) produced similar results.
This random effect structure has been specified and provided to a simr modelling function to simulate the models that will be fitted on our data.
Specifying the random effect structure and the final generative model
## Fixed effects: initialized at 0fixed <-c(700, # intercept0,0,0, # main effects0,0,0, # 2-way interactions0# 3-way interactions )/1000## Random intercepts for participantsv1 <-0.150## Random intercept, slopes (main effects only) and their covariancesslope <- .05correl <-0v2 <-matrix(c(1, correl, correl, correl, 1, correl, correl, correl, 1),3)v2 <-cor_to_cov(v2, sd =c((v1), slope, slope))## residual std devres <- .15model <-makeLmer(formula = y ~ (group + congruence + color)^3+ (congruence + color|id),fixef = fixed, VarCorr =list(v2), sigma = res,data = df )
3 Simulations
Effect sizes: For the interaction Group x Congruence, the effect size can be defined as the double difference (aphantasics - controls) x (congruent - incongruent). Assuming, based on previous works, that the strongest difference would not exceed 50-60ms, and that differences below 10ms would hardly be meaningful, we ran simulations for differences between 10 and 60ms, by steps of 5ms. As the power reached the ceiling of 100% after 40ms effect sizes, we report here the results up to this size.
Number of simulations: For each effect size, we ran as many simulations as were needed to obtain a confidence interval below ±10 points around the mean. Going by steps of 50, we eventually ran 200 simulations for each effect size. These simulations ran for more than 8 hours, therefore their results are stored in .RDS R objects for ease of access: data/r-data-structures/power-analyses-list.RDS and data/r-data-structures/power-analyses-table.RDS. They do not run automatically when rendering the .qmd file.
Running the simulations for various effect sizes
# ─── Running simulations with parallel processing ─────────────────────────────# finding the available cores for parallel processingn_cores <- parallel::detectCores() -1# creating the cluster of coresparallel_cluster <- parallel::makeCluster( n_cores,type ="PSOCK" )# registering the cluster for `foreach`doParallel::registerDoParallel(cl = parallel_cluster)# checking# foreach::getDoParRegistered()# foreach::getDoParWorkers()# initializing, ready for takeoffnsim <-200# Run calculations for a range of effect sizes(simulations <-foreach(effectsize =c(10, 15, 20, 25, 30, 35, 40), .packages =c("tidyverse", "simr") ) %dopar% {# Calculate regression coefficients corresponding to the effect size and # update the model to simulate from accordinglyfixef(model) <-extract_coefficients(effectsize)# Simulate data and estimate powerpowerSim( model, nsim = nsim,test = simr::fixed(xname ="group:congruence", method ='anova' )) }) |>system.time() -> time_simulating# stopping the cluster when we're doneparallel::stopCluster(cl = parallel_cluster)
4 Results
The analysis yielded a power of 82% starting at effect sizes around 25ms, going up to 91% at 30ms and 97% at 35ms. Every effect size equal and above 40ms reached a ceiling of 100% power. The complete results are summarized in Table S4.1 and Figure S4.1.
Generating the summary table for the power analysis
simulations_table <-readRDS(file ="data/r-data-structures/power-analyses-table.RDS")simulations_table |># All the code below is simple formatting of the table that is in the RDSrename("Effect size (difference in ms)"= effectsize,"Successes"= successes,"Simulations"= trials,"Power"= mean ) |>mutate(across(c(lower:upper), ~round(.x, digits =2))) |>unite("CI", lower:upper, sep =", ") |>mutate(CI =paste("[", CI, "]")) |>rename("95% CI"= CI) |>display()
Table S4.1: Results of the power analysis through simulation.
Effect size (difference in ms)
Successes
Simulations
Power
95% CI
10
46
200
0.23
[ 0.17, 0.29 ]
15
78
200
0.39
[ 0.32, 0.46 ]
20
136
200
0.68
[ 0.61, 0.74 ]
25
164
200
0.82
[ 0.76, 0.87 ]
30
182
200
0.91
[ 0.86, 0.95 ]
35
195
200
0.97
[ 0.94, 0.99 ]
40
200
200
1.00
[ 0.98, 1 ]
Generating the plot of the power analysis results
simulations_table |>ggplot(aes(x = effectsize,y = mean,ymin = lower,ymax = upper,label = effectsize )) +# the band representing the 95% CIsgeom_ribbon(fill ="aquamarine", alpha = .2) +# the means and CIs at each specific effect sizegeom_pointrange(colour ="aquamarine4") +geom_point(colour ="aquamarine4") +# line connecting the dotsgeom_line(colour ="aquamarine4") +# horizontal line at 80% power and corresponding annotationgeom_hline(yintercept = .8, color ="black", linetype =1) +annotate(geom ="text",label ="80% power",fontface ="italic",color ="black",x =-Inf,y = .8,vjust =-.5,hjust =-.3,size =4 ) +# horizontal dotted line at 90% power and corresponding annotationgeom_hline(yintercept = .9, color ="gray60", linetype =2) +annotate(geom ="text",label ="90% power",fontface ="italic",color ="grey60",x =-Inf,y = .9,vjust =-.5,hjust =-.3,size =4 ) +# Aestheticsscale_x_continuous(name ="Effect size (RT difference in ms)",breaks =unique(simulations_table$effectsize) ) +scale_y_continuous(name ="Statistical power",breaks =seq(0, 1, .1),limits =c(0, 1) ) +theme_modern() +theme(panel.grid.major =element_line())
Figure S4.1: Results of the power analysis using simulations. 200 simulations have been computed for each effect size.
Figure S4.1: Results of the power analysis using simulations. 200 simulations have been computed for each effect size.
References
Brysbaert, M., & Stevens, M. (n.d.). Power analysis and effect size in mixed effects models: A tutorial. Journal of Cognition, 1(1), 9. https://doi.org/10.5334/joc.10
Fucci, E., Abdoun, O., Baquedano, C., & Lutz, A. (2023). Ready to help, no matter what you did: Responsibility attribution does not influence compassion in expert Buddhist practitioners. Preprint. https://doi.org/10.31234/osf.io/d6y4w